function KF0(H_t,Se,Phi1,y,n_agg)

(T,n) = size(y)
s     = zeros(T,n);
P     = zeros(T,n,n);

# Initialization
s[1,:]   = y[1,:];
P[1,:,:] = H_t[:,:,1];;

# Kalman Filter Recursion
sprime   = zeros(n,1);
liki     = zeros(T,1);

# Forward Iterations
for t=2:T

    # Updating Step
    sprime = Phi1*s[t-1,:];
    Pprime = Phi1*P[t-1,:,:]*Phi1' + Se;

    # prediction step
    yprediction = sprime;
    v = y[t,:] - yprediction;
    F = Pprime + H_t[:,:,t];
    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    liki[t]    = -0.5*n*log(2*pi) - 0.5*log(det(F)) - 0.5*v'*inv(F)*v;

    # updating step
    kgain    = Pprime*inv(F);
    s[t,:]   = sprime + kgain*v;
    P[t,:,:] = Pprime - kgain*Pprime;

end

return s, liki

end

#-----------------------------------------------------------

#-----------------------------------------------------------

function KFnoX(H_t,Se,Phi1,y,n_agg)

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s);
P     = zeros(T,n_s,n_s);

y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition PHI1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
s[1,:]   = s_obs[1,:];
P[1,:,:] = H_t[:,:,1];

# Kalman Filter Recursion
liki     = zeros(T,1);

# Forward Iterations
for t=2:T

    # forecast state and yobs
    yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*s[t-1,:];
    sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*s[t-1,:];

    Pprime_yy = Phi_ys*P[t-1,:,:]*Phi_ys' + Se_yy;
    Pprime_sy = Phi_ss*P[t-1,:,:]*Phi_ys' + Se_sy;
    Pprime_ss = Phi_ss*P[t-1,:,:]*Phi_ss' + Se_ss;

    # condition distr of state on y_obs
    sprime_y_obs    = sprime + Pprime_sy*inv(Pprime_yy)*(y_obs[t,:]-yprime);
    Pprime_ss_y_obs = Pprime_ss - Pprime_sy*inv(Pprime_yy)*Pprime_sy';

    # forecast observation
    yobs_pred = yprime
    sobs_pred = sprime_y_obs;
    v_y       = y_obs[t,:] - yobs_pred;
    v_s       = s_obs[t,:] - sobs_pred;
    F_yy      = Pprime_yy;
    F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    liki[t]  = ( -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                 -0.5*n_s*log(2*pi) - 0.5*log(det(F_ss_yy)) - 0.5*v_s'*inv(F_ss_yy)*v_s)

    # updating step
    kgain      = Pprime_ss_y_obs*inv(F_ss_yy);
    s[t,:]     = sprime_y_obs + kgain*v_s;
    P[t,:,:]   = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

end

return s, liki, P

end


#-----------------------------------------------------------

#-----------------------------------------------------------

function KFwithX(H_t,Se,Phi1,y,n_agg,densdraws,employdraws,knots,xgrid,coef_lambda,coef_mean)

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s);
P     = zeros(T,n_s,n_s);

y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition PHI1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
s[1,:]   = s_obs[1,:];
P[1,:,:] = H_t[:,:,1];

# Kalman Filter Recursion
lh_ss     = zeros(T,1);
lh_X      = zeros(T,1);

# Forward Iterations
for t=2:T

    # forecast state and yobs
    yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*s[t-1,:];
    sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*s[t-1,:];

    Pprime_yy = Phi_ys*P[t-1,:,:]*Phi_ys' + Se_yy;
    Pprime_sy = Phi_ss*P[t-1,:,:]*Phi_ys' + Se_sy;
    Pprime_ss = Phi_ss*P[t-1,:,:]*Phi_ss' + Se_ss;

    # condition distr of state on y_obs
    sprime_y_obs    = sprime + Pprime_sy*inv(Pprime_yy)*(y_obs[t,:]-yprime);
    Pprime_ss_y_obs = Pprime_ss - Pprime_sy*inv(Pprime_yy)*Pprime_sy';

    # forecast observation
    yobs_pred = yprime
    sobs_pred = sprime_y_obs;
    v_y       = y_obs[t,:] - yobs_pred;
    v_s       = s_obs[t,:] - sobs_pred;
    F_yy      = Pprime_yy;
    F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    lh_ss[t]  = ( -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                 -0.5*n_s*log(2*pi) - 0.5*log(det(F_ss_yy)) - 0.5*v_s'*inv(F_ss_yy)*v_s)

    # reconstruct likelihood function for cross-sectional observations
    densdraws_t     = densdraws[1:N,t]
    employdraws_t   = employdraws[1:N,t]
    selecteddraws_t = densdraws_t[employdraws_t.==1]
    coef_t          = coefRecover(s_obs[t,:]', coef_lambda, coef_mean)
    lh_X[t]         = -length(selecteddraws_t)*logspline_obj(selecteddraws_t, coef_t', knots, minimum(xgrid), maximum(xgrid))

    # updating step
    kgain      = Pprime_ss_y_obs*inv(F_ss_yy);
    s[t,:]     = sprime_y_obs + kgain*v_s;
    P[t,:,:]   = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

end

return s, lh_ss, lh_X

end


#-----------------------------------------------------------

function KFnoX_MoreLags(H_t,Se,Psi1,y,n_agg,nlags)

# companion form observations start at t=nlags y[1] and run to t=Tstar y[end]
# H_t starts at t=1 because we need the lags to initialize P_ss for t=0
# Here TmT0p1 is T-T0+1=Tstar-nlags+1

(TmT0p1,n) = size(y)
p       = nlags
n_z     = n_agg*p
n_s     = Int(n-n_z)
n_cross = Int(n_s/p)
s       = zeros(TmT0p1,n_s);
P_ss    = zeros(TmT0p1,n_s,n_s);

z_obs = y[:,1:n_z]
s_obs = y[:,n_z+1:end]

# partition Psi1
Psi_zz = Psi1[1:n_z,1:n_z]
Psi_zs = Psi1[1:n_z,n_z+1:end]
Psi_sz = Psi1[n_z+1:end,1:n_z]
Psi_ss = Psi1[n_z+1:end,n_z+1:end]

Se_yy = Se[1:n_agg,1:n_agg]
Se_ya = Se[1:n_agg,n_agg+1:end]
Se_ay = Se[n_agg+1:end,1:n_agg]
Se_aa = Se[n_agg+1:end,n_agg+1:end]

select_y_from_z  = kron([1.0 zeros(1,p-1)], Matrix(1.0I, n_agg,n_agg))
generate_s_from_a = kron([1.0; zeros(p-1,1)], Matrix(1.0I, n_cross,n_cross))

# Initialization
# tt=1 corresponds to period nlags
s[1,:]   = s_obs[1,:];

diag_iniP = 0.0;
for ll=p:-1:1
    diag_iniP = [diag_iniP; diag(H_t[:,:,ll])]
end
P_ss[1,:,:] = diagm(diag_iniP[2:end])

# Kalman Filter Recursion
liki     = zeros(TmT0p1,1);

# Forward Iterations
# tt=2 corresponds to nlags+1,
for t=2:TmT0p1

    # forecast state and yobs
    yprime = select_y_from_z*(Psi_zz*z_obs[t-1,:]+Psi_zs*s[t-1,:]);
    sprime = Psi_sz*z_obs[t-1,:]+Psi_ss*s[t-1,:];

    Pprime_yy = select_y_from_z*(Psi_zs*P_ss[t-1,:,:]*Psi_zs')*select_y_from_z' + Se_yy;
    Pprime_sy = Psi_ss*P_ss[t-1,:,:]*Psi_zs'*select_y_from_z' + generate_s_from_a*Se_ay;
    Pprime_ss = Psi_ss*P_ss[t-1,:,:]*Psi_ss' + generate_s_from_a*Se_aa*generate_s_from_a';

    # condition distr of state on y_obs
    sprime_y_obs    = sprime + Pprime_sy*inv(Pprime_yy)*(select_y_from_z*z_obs[t,:]-yprime);
    Pprime_ss_y_obs = Pprime_ss - Pprime_sy*inv(Pprime_yy)*Pprime_sy';

    # forecast errors
    v_y       = select_y_from_z*z_obs[t,:] - yprime;
    v_a       = generate_s_from_a'*(s_obs[t,:] - sprime_y_obs);
    F_yy      = Pprime_yy;
    F_aa_yy   = generate_s_from_a'*Pprime_ss_y_obs*generate_s_from_a + H_t[:,:,nlags+t-1]

    liki[t]  = ( -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                 -0.5*n_s*log(2*pi) - 0.5*log(det(F_aa_yy)) - 0.5*v_a'*inv(F_aa_yy)*v_a)

    # updating step
    kgain      = Pprime_ss_y_obs*generate_s_from_a*inv(F_aa_yy);

    s[t,:]       = sprime_y_obs + kgain*v_a;
    P_ss[t,:,:]  = Pprime_ss_y_obs - kgain*generate_s_from_a'*Pprime_ss_y_obs;

end

return s, liki, P_ss

end

function KSmoothnoX_MoreLags(s_filtered,P_filtered,Se,Psi1,y,n_agg,nlags)

    (TmT0p1,n) = size(y)
    n_z     = Int(n_agg*p)
    n_s     = Int(n-n_z)
    n_cross = Int(n_s/p)
    z_obs   = y[:,1:n_z]
    s_obs   = y[:,n_z+1:end]

    a_sim = zeros(TmT0p1,n_cross)
    s_sim = zeros(TmT0p1,n_s)

    # matrices that will be used below
    M = zeros(n, Int(n/p))
    M[1:n_agg,1:n_agg] = Matrix(1.0I,n_agg,n_agg)
    M[n_agg*p+1:n_agg*p+n_cross, n_agg+1:n_agg+n_cross] = Matrix(1.0I,n_cross,n_cross)
    Psi1_p    = Psi1^p
    Sig_sum_p = zeros(n,n)
    for hh = 0:(p-1)
        Sig_sum_p = Sig_sum_p + (Psi1^hh)*M*Se*M'*(Psi1^hh)'
    end
    Xi_s = [zeros(n_z,n_s); Matrix(1.0I, n_s, n_s)]

    # Periods T*, T*-1, ..., T*-p+1
    s_sim[TmT0p1,:]            = s_filtered[end,:] + rand(MvNormal(zeros(n_s), Symmetric(P_filtered[end,:,:])))
    a_sim[TmT0p1-p+1:TmT0p1,:] = reverse(reshape(s_sim[TmT0p1,:],n_cross,p)',dims=1)

    # backward recursion in blocks of length p
    n_blocks = TmT0p1 ÷ p
    T_stop   = TmT0p1 % p

    for tt = (TmT0p1-p):-p:1
        P_tp_t = Psi1_p*Xi_s*P_filtered[tt,:,:]*Xi_s'*Psi1_p' + Sig_sum_p;
        w_tp_t = Psi1_p*[z_obs[tt,:]; s_filtered[tt,:]]
        v_w    = [z_obs[tt+p,:];s_sim[tt+p,:]]-w_tp_t

        s_t_tp = s_filtered[tt,:] + P_filtered[tt,:,:]*Xi_s'*Psi1_p'*inv(P_tp_t)*v_w
        P_t_tp = P_filtered[tt,:,:] - P_filtered[tt,:,:]*Xi_s'*Psi1_p'*inv(P_tp_t)*Psi1_p*Xi_s*P_filtered[tt,:,:]

        s_sim[tt,:] = s_t_tp + rand(MvNormal(zeros(n_s), Symmetric(P_t_tp)))
        if tt-p+1 > 0
           a_sim[tt-p+1:tt,:] = reverse(reshape(s_sim[tt,:],n_cross,p)',dims=1)
        else
           a_sim[1:tt,:] = reverse(reshape(s_sim[tt,1:n_cross*tt],n_cross,tt)',dims=1)
        end
    end

    # to deal with initialization re-run for tt=1
    if (T_stop != 1) & (p > 1)
        P_tp_t = Psi1_p*Xi_s*P_filtered[1,:,:]*Xi_s'*Psi1_p' + Sig_sum_p;
        w_tp_t = Psi1_p*[z_obs[1,:]; s_filtered[1,:]]
        v_w    = [z_obs[1+p,:];s_sim[1+p,:]]-w_tp_t

        s_t_tp = s_filtered[1,:] + P_filtered[1,:,:]*Xi_s'*Psi1_p'*inv(P_tp_t)*v_w
        P_t_tp = P_filtered[1,:,:] - P_filtered[1,:,:]*Xi_s'*Psi1_p'*inv(P_tp_t)*Psi1_p*Xi_s*P_filtered[1,:,:]

        s_sim[1,:] = s_t_tp + rand(MvNormal(zeros(n_s), Symmetric(P_t_tp)))
        a_sim[1,:] = s_sim[1,1:n_cross]
    end

    # Splice the pieces together, removing the duplicate value for t=nlags
    a_sim_initial_lags = reverse(reshape(s_sim[1,:],n_cross,p)',dims=1)
    a_sim_adj = [a_sim_initial_lags; a_sim[2:end,:]]

return a_sim_adj, a_sim, s_sim

end
